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Abstract. The effects of viscoelasticity on the dynamics and break-up of fluid threads in microfluidic T-junctions 
are investigated using numerical simulations of dilute polymer solutions at changing the Capillary number (Ca), i.e. at 
changing the balance between the viscous forces and the surface tension at the interface, up to Ca « 3 x 10 -2 . A Navier- 
Stokes (NS) description of the solvent based on the lattice Boltzmann models (LBM) is here coupled to constitutive 
equations for finite extensible non-linear elastic dumbbells with the closure proposed by Peterlin (FENE-P model). 

We present the results of three-dimensional simulations in a range of Ca which is broad enough to characterize all the 
three characteristic mechanisms of breakup in the confined T-junction, i.e. squeezing , dripping and jetting regimes. The 
various model parameters of the FENE-P constitutive equations, including the polymer relaxation time Zp and the finite 
extensibility parameter L 2 , are changed to provide quantitative details on how the dynamics and break-up properties 
are affected by viscoelasticity. We will analyze cases with Droplet Viscoelasticity (DV), where viscoelastic properties 
are confined in the dispersed (d) phase, as well as cases with Matrix Viscoelasticity (MV), where viscoelastic properties 
are confined in the continuous (c) phase. Moderate flow-rate ratios Q « <^(1) of the two phases are considered in 
the present study. Overall, we find that the effects are more pronounced in the case with MV, as the flow driving the 
break-up process upstream of the emerging thread can be sensibly perturbed by the polymer stresses. 

PACS. 47.50.Cd Non-Newtonian fluid flows Modeling - 47.11.St Multi-scale methods - 87.19.rh Fluid transport 
and rheology - 83.60.Rs Shear rate-dependent structure 


1 Introduction 

Droplet-based microfluidic devices have gained a considerable 
deal of attention, due to their importance in studies that require 
control over droplet size d[33m00[71[H. Common droplet 
generator designs used in these devices are T-shaped mm and 
flow-focusing IUllll2lfl3ll geometries. In T-shaped geometries, 
a dispersed (d) phase is injected perpendicularly into the main 
channel containing a continuous (c) phase. Forces are created 
by the cross-flowing continuous phase which periodically pro¬ 
duces break-up of droplets. The operational regime of these de¬ 
vices is primarily characterized by the Capillary number, which 
quantifies the importance of the viscous forces with respect to 
the surface tension forces at the non-ideal interface, and the 
droplet size and its distribution are dictated by the flow-rate ra¬ 
tio Q = Qd/Qc of the two immiscible fluids. Distinct regimes 
of formation of droplets have been identified: squeezing , drip¬ 
ping and jetting , providing a unifying picture of emulsification 
processes typical of microfluidic systems l9ll 10111 211 l4l . The 
squeezing mechanism of break-up is peculiar of all microflu¬ 
idic systems, because of the physical confinement which natu¬ 
rally accompanies these geometries. In this regime, the break¬ 
up process is driven by the build-up of pressure upstream of 


the emerging thread. The dripping regime, while apparently 
homologous to the unbounded case, is also significantly in¬ 
fluenced by the constrained geometry (9), which modifies the 
scaling law for the size of the droplets derived from the balance 
of interfacial and viscous stresses. Finally, the jetting regime 
sets in only at very high flow rates, or with low interfacial ten¬ 
sion, i.e. higher values of the Capillary number. 

With few exceptions rBinMm previous research has been 
mainly restricted to Newtonian fluids. However, the process¬ 
ing of biological fluids inevitably results in considering a non- 
Newtonian viscoelastic behaviour. Consistently, the study of 
viscoelastic liquids in flow-focusing geometries [15EME) or 
T-junction geometries m has gained some attention. The for¬ 
mation and the pinch-off mechanism of viscoelastic droplets in 
Newtonian continuous phases was investigated in various flow- 
focusing geometries by Steinhaus et al. US), while the effect of 
polymer molecular weight on filament thinning was studied by 
Arratia et al. IfTMTsl . In a recent paper, Derzsi et al. 03 pre¬ 
sented an experimental study of the effects of viscoelasticity 
in microfluidic flow-focusing geometries. The authors find that 
the viscoelasticity of the focusing liquid stabilizes the jets fa¬ 
cilitating formation of smaller droplets, and leads to transitions 
between various regimes at lower ratios of flow and at lower 
values of the Capillary numbers in comparison to the Newto¬ 
nian focusing liquids. Complementing these results with sys¬ 
tematic investigations by varying deformation rates and non- 
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Newtonian constitutive parameters would be of extreme inter¬ 
est. This is witnessed by the various papers in the literature (9J 
[10111 llfl2ll 1911201121112211231 1 addressing these kind of problems 
with the help of numerical simulations. 

Here we present a three-dimensional numerical investigations 
of the interplay between viscoelasticity and geometry-mediated 
breaking in confined microfluidic T-junctions. Numerical sim¬ 
ulations allow to address systematically the importance of the 
various free parameters in the viscoelastic model and visualize 
the distribution of the polymer feedback stresses, thus correlat¬ 
ing the distribution of those stresses to the interface shape. Our 
numerical approach offers the possibility to tune the viscosity 
ratio of the two Newtonian phases, a fact that is instrumental to 
perform simulations with non-Newtonian phases and compare 
them with the results of fully Newtonian systems with the same 
viscosity ratio. 

The paper is organized as follows: in Sec. [2] we will present 
the necessary mathematical background for the problem stud¬ 
ied, showing the relevant equations that we integrate in both 
the continuous and dispersed phases, and identifying the rele¬ 
vant dimensionless numbers useful for our investigation. Use¬ 
ful benchmarks for the shear rheology of the numerical model 
will be provided for the typical parameters used in our study. 
In section [3] we will present the numerical results and charac¬ 
terize the effects of viscoelasticity in the three distinct regimes 
of squeezing (subsection HU and dripping/jetting (subsection 


3.2). We will study both the droplet size soon after break-up as 


well as the characteristic time for break-up and compare them 
with the corresponding Newtonian cases. To explain the ob¬ 
served behaviour we will explore the distribution of feedback 
stresses in the non-Newtonian phases and correlate them with 
the characteristic mechanisms of break-up in the confined T- 
junction. Conclusions will follow in section [4] 


2 Theoretical Model 

Numerical modeling of viscoelastic fluids often relies on the 
coupling of constitutive relations for the stress tensor, typically 
obtained via approximate representations of some underlying 
micro-mechanical model for the polymer molecules, with a 
Navier-Stokes (NS) description for the solvent. The FENE-P 
constitutive model is obtained via a pre-averaging approxima¬ 
tion applied to a suspension of non interacting finitely exten¬ 
sible non-linear elastic (FENE) dumbbells. FENE-P is well- 
adapted for dilute (and semi-dilute) polymer solutions, and has 
been used previously to analyze filament thinning of viscoelas¬ 
tic fluids in macroscopic experiments II241I25L as well as the 
effects of viscoelasticity on the dynamics of filament thinning 
and break-up processes in microchannels dm. In this paper 
we provide quantitative details on how the FENE-P model pa¬ 
rameters affect the break-up properties of confined threads in 
microfluidic T-junctions, by analyzing separately the cases of 
Droplet Viscoelasticity (DV), where the viscoelastic properties 
are confined in the dispersed (d) phase undergoing the break-up 
process, as well as the cases with Matrix Viscoelasticity (MV), 
where the viscoelastic properties are confined in the continuous 
(c) phase. A fluid described by the FENE-P model possesses 
the same dynamical properties as a fluid described by the much 
simpler Oldroyd-B model, which assumes that polymers can 


be modeled as Hookean springs which relax to the equilibrium 
configuration with a characteristic time T p. The main differ¬ 
ence is that the Oldroyd-B model allows for infinite extension 
of polymer molecules, while the FENE-P model uses a spring- 
force law in which the polymer molecules can be stretched only 
by a finite amount in the flow field 1126112711 . Thus we can ex¬ 
plore systematically both the effects of the polymer relaxation 
times as well as their finite extensibility. 

The solvent part of the model is obtained with lattice Boltz¬ 
mann models (LBM) [28,29], which proved to be extremely 
valuable tools for the simulation of droplet deformation prob¬ 
lems 130113 lll32ll33lL droplets dynamics in open 13411351 and 
confined [22l|23j33i microfluidic geometries. LBM is instru¬ 
mental to solve the diffuse-interface hydrodynamic equations 
of a binary mixture of two components 136U37U381I39U401I41M : 
the resulting physical domain can be partitioned into different 
subdomains, each occupied by a “pure” fluid, with the inter¬ 
face between the two fluids described as a thin layer where 
the fluid properties change smoothly. The FENE-P constitutive 
equations are solved with a finite difference scheme which is 
coupled with the solvent LBM as described in EUa.Thenu- 
merical approach has been extensively validated in our previ¬ 
ous works mm, where we have provided evidence that the 
model is able to capture quantitatively rheological properties 
of dilute suspensions as well as deformation and orientation of 
single viscoelastic droplets in confined shear flows. The main 
essential features of the model are recalled in Appendix [A| 

In the MV case, the equations we solve in the continuous phase 
are the Navier-Stokes (NS) equations coupled to the FENE-P 
constitutive equations 


Pc \dtUc V (Uc ’ VK] — V Pc~\~ 

V (t} c (Vm c + ('Vm c ) 7 )) + ^ V • [f{rp)C], { 1} 

'T P 

d t C + (u c ■ V)C = C ■ {Vu c )+(Vu c ) T ■ C 

_ ( Jf(r P )C-I ^ (2) 

Here, u c and r\ c are the velocity and the dynamic viscosity 
of the continuous phase, respectively. p c is the solvent den¬ 
sity, P c the solvent bulk pressure, and (Vw c ) r the transpose 
of (Vit c ). As for the polymer details, r)p is the viscosity 
parameter for the FENE-P solute, Tp the polymer relaxation 
time, C the polymer-conformation tensor, I the identity tensor, 
f{rp) = (L 2 — 3)/(L 2 — rp) the FENE-P potential that ensures 
finite extensibility, rp = yjTr(C) and L is the maximum pos¬ 
sible extension of the polymers 12611271 . In the dispersed phase 
we just consider the NS equations 

Pd [d t u d + (u d • V)u d ] = -VP d 

+V(i)j(V% + (V«/)) (3) 

where the different fields have the same physical meaning but 
they refer to the dispersed phase. Immiscibility between the 
dispersed phase and the continuous phase is introduced using 
the so-called “Shan-Chen” model Il42ll44ll45l which ensures 
phase separation with the formation of stable interfaces be¬ 
tween the two phases characterized by a positive surface ten¬ 
sion a. 
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For the DV case, we consider the reversed case, where the 
FENE-P constitutive equations are integrated in the dispersed 
phase (i.e. 0-0 with c —y d), while only the NS equations are 
considered in the continuous phase (i.e. 0 with d -A c). 

As for the geometry used, the T-junction is embedded in a rect¬ 
angular parallelepiped with size L x xL y xL z , and channels have 
a square cross-section with edge H = L z . The square cross- 
section is resolved with H xH = 32x32 grid points. The main 
channel and the side channel lengths are resolved with a vari¬ 
able number of grid points (see also table [TJ, depending on the 
characteristic regime analyzed and the characteristic size of the 
droplet after break-up. 

Besides the geometrical parameters, the Newtonian problem is 
described by six parameters characterizing the flow and ma¬ 
terial properties of the fluids. These parameters are the mean 
speeds of the continuous and dispersed phases, v c and vj, re¬ 
spectively; the viscosities of the two fluids r\ c and of Eqs. 
0 and 0, the interfacial tension cr, and the total density p c = 
pd = p (the same for the dispersed and continuous phases). We 
will assume perfect wetting for the continuous phase, while the 
dispersed fluid does not wet the walls. Wetting properties are 
introduced at the boundaries declaring the stress of the density 
fields B5B7I . We then choose the following groups mm2: 
the Capillary number calculated for the continuous phase, 


Ca = 



(4) 


the Reynolds number Re = pv c ///(r] TOT?c ), the viscosity ratio 
A, and the flow rate ratio 


Q = 


Vd 

Vc 


Qd 

Qc 


(5) 


where Qd = VdH 2 and Q c = v c H 2 are the flow rates at the two 
inlets. For the flow regimes under consideration, the Reynolds 
number is small (Re ~ 0.01 — 0.1), and does not influence the 
droplet size, which leaves us with the three governing param¬ 
eters: Ca, A and Q. Notice that the total viscosity in the con¬ 
tinuous phase T7 tot?c is either 7] TOT?c = flc + (for MV) or 
Atot,c — flc (for DV). In the outlet, we impose pressure bound¬ 
ary conditions and use Neumann boundary conditions for the 
velocity field. A Dirichlet boundary condition is imposed at 
the inlets by specifying the pressure gradient that is compat¬ 
ible with the analytical solution of a Stokes flow in a square 
duct ED As for the polymer boundary conditions, we impose 
a Dirichlet type boundary conditions by linearly extrapolating 
the conformation tensor at the boundaries. 

Our numerical approach offers the possibility to tune the vis¬ 
cosity ratio of the two Newtonian phases H421I43I1 . This will al¬ 
low us to work with unitary viscosity ratio, defined in terms of 
the total (fluid + polymer) shear viscosity A = + fl/ 5 ) = 

1.0 for MV and A = ( rjd + r lp)/ r lc = 1.0 for DV. Consis¬ 
tently, we will compare the non-Newtonian simulations with 
the corresponding Newtonian case at A = = 10. The 

ratio between the polymer viscosity and the total viscosity is 
set to rip/(r[ Ci d + Vp) ~ 0.265. Similarly to problems involv¬ 
ing single droplet deformation and dynamics 148,49,50"51 
152), we choose to quantify the degree of viscoelasticity with the 


Deborah number that we define as De = 


NjH ( o \ 
\{r\ diC +r]p)H1) ’ 


where N\ is the first normal stress difference which develops 
in the viscoelastic phase in presence of a homogeneous steady 
shear 12611251 . In the definition of the Deborah number, the vis¬ 
cosity is obviously indicated in the viscoelastic phase, either 
r\ c + Vp f° r MV or rid+ rip for DV. The shear rheology of the 
model can be quantitatively verified in the numerical simula¬ 
tions. There are indeed exact analytical results one can get by 
solving the constitutive equations for the hydrodynamical prob¬ 
lem of steady shear flow, u x = yy, u y = u z = 0: both the poly¬ 
mer shear stress and the first normal stress difference N\ for the 
FENE-P model 1261125) follow 


Tf('-rKy = 

Tp 


2ljp 

Tp 



sinh 


^arcsinh 




( 6 ) 


Ni 


rip (L 


^/(/>)(<:„ Cyy) - 8 ^ g 


sinh 2 


-arcsinh 

3 




(7) 


The validity of both Eqs. ([6]) and ([7]) is benchmarked in Fig. 
[TJ numerical simulations have been carried out in three dimen¬ 
sional domains with H x H x H = 20 x 20 x 20 cells. Periodic 
boundary conditions are applied in the stream-flow (x) and in 
the transverse-flow (z) directions while two walls are located 
at y = 0 and y = H. The linear shear flow u x = yy, u y = u z = 0 
is imposed in the numerics by applying two opposite veloci¬ 
ties in the stream-flow direction ( u x (x,y = 0,z) = —u x (x,y = 
H.z) = U w ) at the upper (y = H) and lower wall (y = 0) with 
the bounce-back rule E2 We next change the shear in the 
range 10 -6 < 2 U w /H < 10 -2 lbu (lattice Boltzmann units) and 
the polymer relaxation time in the range 10 1 < Tp < 10 5 lbu 
for different values of the finite extensibility parameter ranging 
in the interval L 2 = 5 x 10 — 5xl0 3 , and fixed T]p. The var¬ 
ious quantities are made dimensionless with the viscosity T]p 
and the relaxation time Tp, and they are plotted as a function 
of the dimensionless shear A = Tpj. The values of the con¬ 
formation tensor are taken when the simulation has reached 
the steady state. As we can see from the figures, all the nu¬ 
merical simulations collapse on different master curves, depen- 
dently on the value of L 2 . In particular, both the stress 0 and 
first normal stress difference ([7]) increase at large A to exhibit 
variable levels depending on Lr, and consistently with the the¬ 
oretical predictions H261I271I251 . The dependence from L 2 re¬ 
flects in thinning effects visible in the dimensionless polymer 
shear viscosity, f(rp)C xy /A, and first normal stress coefficient, 
= f(rp)(C xx — C yy )/A 2 , which are analyzed in the bottom 
panel of figure[l] Overall, the numerical simulations performed 
to quantify the shear rheology reveal a very good agreement 
with the theoretical predictions both in the polymer shear vis¬ 
cosity and in the first normal stress difference. Similar analysis 
can be performed for extensional flows, showing that the in¬ 
crease of the extensional viscosity predicted by the theory [26, 
1271251 is indeed found in the numerical simulations 142) . The 
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coupling between normal stresses and single droplet dynam¬ 
ics under simple shear has also been extensively verified in the 
numerical simulations. In particular, in lf42l we provided ev¬ 
idence that the model proposed captures quantitatively single 
droplet orientation and deformation in presence of viscoelastic 
stresses. 

In the limit of Hookean dumbbells (Oldroyd-B limit, L 2 » 1) 
we can use the asymptotic expansion of the hyperbolic func¬ 
tions and we get N\ « 2zprjpy 2 , so that 


De = 


r P r\p 

r\d, c + Ip ’ 


( 8 ) 


Equation ([8]) shows that De is clearly dependent on the ratio 
between the polymer relaxation time Tp and the time t h defined 
as 

„ H (r] dtC + tip) 

= - (9) 

G 

which represents the relaxation time of a droplet with char¬ 
acteristic size H , determined by viscous and capillary forces. 
Clearly, definition is dependent on rheology and geometry. 
The values of L 2 we use in the numerical simulations of the 
confined T-Junctions are such that L 2 > 10 2 , ruling out impor¬ 
tant thinning effects for the shears achieved in our simulations. 
We therefore choose to report results based on the definition 
of the Deborah number ([8]) together with the finite extensibil¬ 
ity parameter L 2 . All the various parameters are summarized 
in Table [I] An interesting point of discussion emerges from 
the attempt of connecting results from numerical simulations 
with experimental data, and in particular how appropriate is 
the choice of the parameters rfp, t p and L 2 . Some of these in¬ 
formation are available from the literature (see 11511181 ! and ref¬ 
erences therein). Arratia et al. EE!) performed experiments 
on filament thinning and break-up of viscoelastic fluids in mi¬ 
crochannels: for a viscoelastic fluid made by adding 100 ppm 
of polyacrylamide (PAA) with MW (molecular weight) of 10 5 , 
a concentration of r\pl[v\d^c + Vp) ~ 10 -1 , a finite extensibil¬ 
ity parameter L 2 « 10 3 and fluid relaxation time Tp = 0.05 s are 
found to best fit the experimental rheological data. The polymer 
relaxation time decreases at decreasing the molecular weight, 
down to t p « lO -3 ^, for MW of 1 x 10 3 . In the present study, 
we choose to use different L 2 , so as to study the enhancement 
of viscoelastic effects up to the value above cited. As for the 
polymer relaxation time Tp, we notice that a Newtonian droplet 
with characteristic size of the order of 10 -4 m would result in 
a T# (Vd ~ 0.2 Pas and a = \0~ 2 N/m 1 15118 0 of the order 
of T h = r\dH/o « 10 -3 v, hence t p /t h ranges from 1 to a few 
tens. Such a range can actually be explored in the numerics 
by tuning T p in the range 250 — 4000 lbu (Tp/t# in the range 
1-25). 


3 Results and Discussions 

In figure [2] we report 3D snapshots illustrating geometry me¬ 
diated break-up in various scenarios depending on Ca. These 
snapshots allow us to identify the various regimes which are 
known from the literature on droplet formation in confined T- 
junctions (see l9lll 1111211 and references therein): these will be 


used as “reference” Newtonian scenarios to quantify the im¬ 
portance of viscoelasticity. Notice that we have used the char¬ 
acteristic shear time T shear = H/v c as a unit of time. At low Ca 
(Panels (a)-(d) of figure|2]), the incoming thread tends to occupy 
and obstruct the entire cross-section of the main channel, with 
the break-up occurring at the junction (Panel (d) in figure [2]). 
By increasing Ca, a dripping scenario is entered (Panels (e)- 
(h) in figure [2]) in which the obstruction of the cross-section 
in the main channel is less visible and viscous shear forces 
start to influence the droplet break-up process immediately af¬ 
ter the droplet enters into the main channel (Panel (f) in figure 
As a result of the combined effect of surface tension and 
viscous forces, smaller droplets are formed downstream of the 
T-junction (see Panels (g)-(h) in figure|2]). By further increasing 
Ca, a critical value 0 exists above which the dispersed phase 
develops a thread entering the main channel and the droplet 
detachment point gradually moves downstream, until a jet is 
formed. The length of the jet is obviously limited by the size 
of the computational domain and simulations with large res¬ 
olution are indeed necessary (see table [I]) to make sure that 
the finite simulation domain does not play a role in the droplet 
formation inside the junction. A quantitative analysis on the in¬ 
fluence of viscosity ratio and channel geometries on the above 
described physical scenarios has already been provided in the 
literature iTOirnt . Here, instead, we aim to illustrate the ef¬ 
fects of viscoelasticity. As already stressed in section [2] our 
numerical approach offers the possibility to tune the viscos¬ 
ity ratio of the two Newtonian phases 1421 431. By fixing the 
polymer viscosity T]p, we can use such flexibility to achieve 
unitary viscosity ratio, defined in terms of the total (fluid + 
polymer) shear viscosity A = r\d/(ri c + t]p) = 1.0 for MV and 
A = (rid + rip) jr( c = 1.0 for DV. This will allow us to compare 
the non-Newtonian simulations with the corresponding Newto¬ 
nian case at the same (unitary) viscosity ratio. We will explore 
systematically both the effects of the finite extensibility param¬ 
eter L 2 and the polymer relaxation time Tp. 


3.1 Squeezing Regime 

To go deeper and be more quantitative on the characterization 
of the various regimes, we start by investigating the droplet size 
as a function of the flow-rate ratio Q in the squeezing regime. 
The characteristic droplet size Lj in the squeezing regime is 
only weekly affected by the viscosity ratio and mainly deter¬ 
mined by the ratio of the volumetric flow-rates of the two im¬ 
miscible fluids as 


Qd 

Ld — ofl + tx 2 ~r~ — cc\-\- a 2 Q. 
He 


( 10 ) 


The constants oq and a 2 , which are of the order one, are de¬ 
termined by the channel geometry lfl2ll . The linear scaling 
law (T0| ) has already been verified in experiments I3ll6lfl4ll 
and also in numerical simulations ll9ll 11 II 12115411 . Our Newto¬ 
nian data in the squeezing regime are quantitatively analyzed 
in figure [3] where we report the dimensionless droplet volume 
V/H 3 = Ld/H. The linear behaviour of Eq. (10) is indeed re¬ 
produced by our simulations (oq = 1 and a 2 = 2) which are 
well in agreement with other existing numerical data in the 
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(a) Polymer shear stress 


(b) Polymer first normal stress difference 




A A 


(c) Polymer shear viscosity 


(d) Polymer first normal stress coefficient 


Fig. 1: Polymer shear rheology. Panels (a)-(b): we plot the polymer shear stress and the first normal stress difference (both scaled 
to the polymer viscosity 7]p and polymer relaxation time Tp, see and text for details) as a function of the dimensionless 

shear A = T>y in a steady shear flow with intensity y. Symbols are the results of the numerical simulations H42U43I1 with different 
imposed shears, different Tp and different L 2 . All the numerical results collapse on different master curves, dependently on the 
value of L 2 : L 2 = 5 x 10 (squares), L 2 = 5 x 10 2 (circles), L 2 = 5 x 10 3 (triangles). The lines are the theoretical predictions 
based on Eqs. ([6]) and ([7]). Panels (c)-(d): we plot the dimensionless polymer shear viscosity and the first normal stress coefficient 
extracted from data in the top panels. 


Ca 

Q 

L x xL y x L z 
cells 

Vd 

lbu 

Ac 

lbu 

Vp 

lbu 

Tp 

lbu 

De 

Lr 

0.002-0.02 

1.0 

640 x 128 x 32 

0.49 

0.49 

0.00 




0.002-0.02 

1.0 

896 x 128 x 32 

0.36 

0.49 

0.13 

2-45 x 10 2 

0.3-7.0 

10 2 ,10 3 ,10 4 

0.002-0.02 

1.0 

896 x 128 x 32 

0.49 

0.36 

0.13 

2-45 x 10 2 

0.3-7.0 

10 2 ,10 3 ,10 4 

0.002 

0.25-1.0 

640 x 128 x 32 

0.49 

0.49 

0.00 




0.002 

0.25-1.0 

896 x 128 x 32 

0.36 

0.49 

0.13 

2-45 xlO 2 

0.3-7.0 

10 2 ,10 3 ,10 4 

0.002 

0.25-1.0 

896 x 128 x 32 

0.49 

0.36 

0.13 

2-45 x 10 2 

0.3-7.0 

10 2 ,10 3 ,10 4 


Table 1: Parameters for the numerical simulations: Ca is the Capillary number (see Eq. ([4j), Q = Qd/Qc is the flow-rate ratio between the 
dispersed (d) and continuous (c) phase. The T-junction is embedded in a rectangular parallelepiped with size L x x L y x L z , and channels have 
a square cross-section with edge H = L z . rfd is the dynamic viscosity of the Newtonian solvent inside the dispersed phase, r/ c is the dynamic 
viscosity of the Newtonian solvent inside the continuous phase, pp is the polymer viscosity, Tp is the polymer relaxation time, De is Deborah 
number based on definition {8}. 
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to + 5.6T sh ear, Q = 0.5, Ca = 0.0026 



(e) t = to + 2.6T shear ,<2= 1.0, Ca = 0.013 




?o + 7.4T shear , Q = 0.5, Ca = 0.0026 



(g) t — ? 0 + 4.0T shear , (2 = 1.0, Ca = 0.013 




(i) t = *o + 2.8T S hear, <2 = 1-0, Ca = 0.026 



(j) t = ?o + 3.8T S hear, <2= 1-0, Ca = 0.026 



(k) t = to + 5.2T shear , Q = 1.0, Ca = 0.026 



(1) t = to + 6.6T s hear 5 Q= 1.0, Ca = 0.026 


Fig. 2: Droplet formation in T-junction geometries for a Newtonian case with viscosity ratio A = 1.0. Panels (a)-(d): we illustrate 
the squeezing regime at Ca = 0.0026 and flow-rate ratio Q = 0.5: the fluid thread enters and obstructs the main channel and 
break-up is mainly driven by the pressure build-up upstream of the emerging thread (9). Both the dynamics of break-up and 
the scaling of the sizes of droplets are influenced weakly by viscous forces am Panels (e)-(h) show typical features of the 
dripping regime at Ca = 0.013 and Q = 1.0: the break-up process starts to be influenced by the shear forces, although the thread 
still occupy a significant portion of the main channel. This results in smaller droplets formed downstream of the T-junction. 
Panels (i)-(l) report snapshots from the jetting regime at larger Capillary number, Ca = 0.026, and Q = 1.0: the dispersed phase 
develops a thread entering the main channel and the droplet detachment point gradually moves downstream, until a jet is formed. 
To better highlight the jetting regime, the associated figures display a larger portion of the main channel of the T-junction. In 
all cases we have used the characteristic shear time T shear = H/v c as a unit of time, while to is a reference time (the same for all 
simulations). 


literature, obtained with phase field numerical simulations m 
and LBM simulations [541. Notice that the numerical simula¬ 
tions of Bower & Lee f54l are performed with a viscosity ratio 
A =0.02. Nevertheless, their results agree with the others (in¬ 
cluding ours), which is a distinctive feature of the squeezing 
regime, where the droplet size is greatly affected by Q and lit¬ 
tle effect is expected from a change in the fluid properties (i.e. 
change in the viscosity ratio A). 

To proceed further, we compute the droplet size for the two dis¬ 
tinct cases of MV and DV. Panel (a) in figure [4] refers to a case 
with MV, with flow-rate ratio and finite extensibility parame¬ 
ter ranging in the interval Q = 0.2 — 1.0 and L 2 = 10 2 — 10 3 , 
respectively. For the non-Newtonian cases, the polymer relax¬ 
ation time has been kept fixed to = 4000 lbu: this is a value 
at which the characteristic Deborah number ([8]) is of order 1 
and viscoelastic effects are clearly visible. Panel (b) of figure 
[4] reports the same quantities as Panel (a) for a case with DV. 
The scaling relation ( fTO] ), which is peculiar of the Newtonian 
cases, is a result of continuity. The analysis of the droplet size 
as a function of the flow-rate ratio Q reveals that such relation 


needs to be modified to account for the effects of viscoelas¬ 
ticity: for increasing finite extensibility parameters, the droplet 
size is manifestly decreased by matrix viscoelasticity. Overall, 
figure]?] conveys the message that viscoelastic effects are more 
pronounced in the case of MV, whereas cases with DV only 
show smaller deviation with respect to the Newtonian refer¬ 
ence case. This is not surprising, in view of the fact that the 
break-up process in the squeezing regime is driven by the ac¬ 
tion of the flow upstream of the emerging thread. More quan¬ 
titatively, the linear scaling law © is the result of two dis¬ 
tinct physical processes: first, the dispersed phase grows until 
it effectively blocks the crosssection of the main channel and 
obstructs the flow of the continuous fluid (see also Panel (a) 
in figure [2]). At this particular moment, the “blocking length” 
Lbiock is of the order of the channel width, say CC\ H (with CC\ a 
constant of order unity). Afterwards, the increased pressure in 
the continuous phase begins to squeeze the neck of dispersed 
phase (see also Panels (b)-(d) in figure [2]). For a neck with 
a characteristic width ( OC 2 is a constant, again, of order 
unity) and squeezing at a rate approximately equal to the av- 
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Fig. 3: Panels (a)-(c): Effect of the flow-rate ratio Q in the squeezing regime with Ca = 0.0026 and A = 1.0. In Panel (d) we 
report the dimensionless droplet volume as a function of the flow-rate ratio Q. Our data are compared with the phase field 
numerical simulations of De Menech et al. 0 and the LBM simulations of Bower & Lee (54). Superimposed we report the 
linear fit predicted by Garstecki et al. fl4) (see Eq. (T0|)), based on the assumption that the droplet size is greatly determined by 
the ratio of the volumetric flow-rates of the two immiscible fluids. Notice that the numerical simulations of Bower & Lee ED 
are performed with a viscosity ratio A = 0.02 which differs from ours. However, in the squeezing regime good agreement is still 
found, since the droplet size is greatly affected by Q and little effect is expected from a change in the fluid properties (i.e. change 
in A). To test the robustness of our findings at changing t he ch annel dimensionality, we repeated the numerical simulations in a 
2d channel with viscosity ratio A =0.05 (see also section 3.2 for discussions). 



(a) Matrix Viscoelasticity (MV), Ca = 0.0026 (b) Droplet Viscoelasticity (DV), Ca = 0.0026 

Fig. 4: Quantitative analysis of the break-up process in the squeezing regime at Ca = 0.0026. We report the dimensionless droplet 
volume V/H 3 soon after break-up for a case with matrix viscoelasticity (MV) and droplet viscoelasticity (DV). We choose the 
flow-rate ratio Q and finite extensibility parameter L 2 ranging in the interval Q = 0.2 — 1.0 and L 2 = 10 2 — 10 3 , respectively. For 
the non-Newtonian cases, the polymer relaxation time has been kept fixed to Tp = 4000 lbu, corresponding to a Deborah number 
De = 5.7, based on definition (8}. Data for different Tp at fixed flow-rate ratio Q = 1.0 are reported in figure [5] 


erage velocity ( Q c /H 2 ), it takes a time T squeeze « 0 C 2 HH 2 /Q c 
to complete the squeezing process. During this time, the thread 
continues to elongate at rate Qd/H 2 . The resulting “squeez- 
ing length” is therefore L squeeze «t sqmeZ eQd / H 2 = a 2 HQ d /Q c . 
Consequently, the final dimensionless size L^/H of the droplet 
can be expressed as L^/H « + 0 C 2 Q. Panel (a) of figure[4]ac- 

tually reveals a change in the slope at increasing L 2 : while the 


slope at L 2 = 100 is still almost same to that of the Newtonian 
case, the slope at L 2 = 1000 is visibly different. This points to 
the fact that the largest elastic effects may effectively perturb 
the region of the fluid upstream of the junction. 

To better complement the results of figure [4| in figure [5] we 
study the droplet size for the same values of L 2 analyzed 
in figure [4] and different values of Tp ranging in the interval 
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Fig. 5: Quantitative analysis of the break-up process in the squeezing regime at Ca = 0.0026. We report the dimensionless droplet 
volume V /V Newt soon after break-up for a case with matrix viscoelasticity (MV). The droplet volume has been made dimensionless 
with respect to the Newtonian volume (V Newt ) for the same Ca. The finite extensibility parameter L 2 and the polymer relaxation 
time t p are ranging in the interval L 2 = 10 2 — 10 3 and T p = 250 — 4000 lbu, respectively. Correspondingly, the Deborah number 
is reported. In all cases, the flow-rate ratio Q and the viscosity ratio between the two fluids have been kept fixed to Q = A = 1.0. 


Tp = 250 — 4000 lbu, resulting in a Deborah number ranging in 
the interval De = 0.4 — 6.2. For the all L 2 studied, the droplet 
size shows a decreasing behaviour at increasing the Deborah 
number, which is more pronounced at larger L 2 . Consistently 
with the expectations, when De ^ 0 we observe minor de¬ 
viations with respect to the Newtonian case. We notice that 
the same analysis (data not shown) for DV reveals only a mi¬ 
nor effect of non-Newtonian rheology in the dispersed phase, 
stressing once more the fact that viscoelastic effects in the up¬ 
stream of the emerging thread are more efficient in perturbing 
the break-up process. 

That viscoelastic effects are more pronounced in presence of 
larger L 2 is qualitatively understood because, by increasing L 2 , 
the polymer dumbbell becomes more extensible and the max¬ 
imum level of stress attainable is increased 1261271 . Consis¬ 
tently, we expect an increased effect of the polymer feedback 
stresses on the Newtonian solvent. However, results of figures 
[4][5] only support this statement indirectly, i.e. without any in¬ 
formation on the distribution of polymer feedback stresses and 
their action on the droplet formation process. To go deeper into 
this point, in figure [6] we report a simultaneous view of the 
droplet shape just before break-up and the polymer feedback 
stress that develops in the non-Newtonian phase. In particular, 
we focus on the polymer feedback stress in the stream-flow di¬ 
rection 

Txx = —f(rp)C xx . (11) 

T P 

We observe that T xx is enhanced in the region upstream of 
the emerging thread, providing extra viscoelastic forces which 
combine to change the droplet break-up process. 

To make progress, we have monitored the time evolution of the 
pressure in the continuous fluid immediately upstream of the 
T-junction. In Panels (a)-(b) of figure [7] we report the pressure 
as a function of time for Ca = 0.0026, with L 2 = 100 (Panel 


(a)), L 2 = 1000 (Panel (b)), and different De. It is evident that 
the obstruction of the main channel leads to an increase of the 
pressure upstream of the T-junction. As the dispersed phase en¬ 
ters the junction, the pressure rises gradually until the channel 
is blocked. The presence of viscoelasticity actually proves in¬ 
strumental to enhance the pressure build-up and the effect is 
more pronounced at increasing both De and L 2 . One may at¬ 
tempt to explain the observed behaviour by arguing that the 
obstruction provided by the thread forces the viscoelastic ma¬ 
trix fluid to “converge” and flow into a constriction, hence to 
develop a high extensional viscosity 12611271. This viscous re¬ 
sponse increases the dissipation and hence the pressure drop. 
This interpretation is actually supported by a direct observation 
of the velocity streamlines in the moment of the obstruction as 
reported in figure [8] Moreover, we have analyzed the force bal¬ 
ance in the whole region upstream of the emerging thread. In 
particular, we have defined an effective force (F eff ) as l43l 

F eB = ^ V • [f(r P )C] - V (tjpO Vu c + (Vu c ) T )) . (12) 

T P 

Indeed, we remark that viscoelastic forces provide a contribu¬ 
tion to the shear forces. This happens in simple shear flows 
and also for weak viscoelasticity [26,27], where we expect that 
the viscoelastic stresses closely follow the viscous stresses, i.e. 

■ L f(rp)C ] » V • (77p(Vw c + {Vu c ) t )). Obviously, this 
cannot be the case when viscoelasticity is enhanced and the 
Deborah number is above unity. Since all our simulations are 
performed with the same shear viscosity, the effective force 
gives an idea of how much the viscoelastic system differs from 
the corresponding Newtonian system with the same viscosity. 
If present (F eff 7 ^ 0), this change is attributed to viscoelasticity. 
For a case with L 2 = 5 x 10 3 , we analyze the effective force 
in the xy-plane at z = L z /2 in the moment when the thread ob¬ 
structs the main channel. Results are reported in figure|9] where 
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Fig. 6: Panels (a)-(c): density contours of the dispersed phase overlaid on the polymer feedback stress in the stream-flow direction 
(see Eqs. Q and ( |11| )) for a case with matrix viscoelasticity (MV) with three different L 2 : L 2 = 10 2 (Panel (a)), L 2 = 5 x 10 2 
(Panel (b)) and L 2 = 10 3 (Panel (c)). All the other parameters are kept fixed: De = 5.7, Ca = 0.0026, A = 1.0 and Q = 1.0. As L 2 
is increased, we see that the flow in the continuous phase develops enhanced polymer feedback stresses upstream of the emerging 
thread. In all cases we have used the characteristic shear time T shear = H/v c as a unit of time, while to is a reference time (the same 
for all simulations). Notice that the colorbar of the feedback stress © is the same. 


we show two distinct plots for the x and y component of F eS . 
Upstream of the emerging thread, and close to the bottom wall, 
we indeed observe a resistance force which opposes to the flow¬ 
ing through the constriction, and we believe is responsible for 
the build-up in the pressure. 

The tendency of viscoelastic stresses to promote a smaller 
droplet volume soon after break-up may be provisionally 
thought of as an anticipation of the dripping regime, thus echo¬ 
ing the work by Derzsi et al. ED in the flow-focusing geom¬ 
etry, where the authors found that the viscoelasticity leads to 
transitions between various regimes at lower ratios of flow and 
at lower values of the Capillary numbers in comparison to the 
Newtonian focusing liquids. However, upon entering the drip¬ 
ping regime, viscous shear forces will become relevant and 
since the shear viscosity is kept the same in all the simulations, 


one should expect to find a less pronounced effect of viscoelas¬ 
ticity at larger Capillary numbers. These expectations are in¬ 
deed borne out by numerical simulations in the next section. 


3.2 Dripping and Jetting Regimes 

The analysis in the squeezing regime has evidenced the non 
trivial role of the polymer feedback stresses in changing the dy¬ 
namics and break-up properties in a situation where Ca is mod¬ 
erately small. Consequently, an interesting point of discussion 
emerges on the role of viscoelasticity on scenarios which are 
different from the squeezing regime. As we have seen in figure 
[2] by increasing Ca at fixed flow-rate ratio we move from the 
squeezing regime to the dripping and jetting regimes. Also, as 
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(a) L 2 = 100, Ca = 0.0026 


(b) L 2 = 1000, Ca = 0.0026 


Fig. 7: Analysis of the pressure (P) versus time in the squeezing regime. The pressure is computed upstream of the T-junction 
and Po is a constant reference pressure computed in the static case. Panel (a): we report the normalized pressure versus time with 
fixed Ca = 0.0026 for the Newtonian case (black squares) and two cases with matrix viscoelasticity (MV) at fixed L 2 = 100: 
De = 1.43 (red circles) and De = 5.7 (blue triangles). Panel (b): same as panel (a) with L 2 = 10 3 . In all cases we have used the 
characteristic shear time T shear =H/v c as a unit of time, while t ohs is the time when the thread starts to obstruct the channel. 



Fig. 8: Left Panel: Velocity streamlines (black lines) overlaid on the polymer feedback stress in the streamflow direction (see Eqs. 
0 and fTT] )) for a case with matrix viscoelasticity (MV). The obstruction provided by the thread forces the flow to converge into 
the gap, hence triggering an extensional response in the fluid region upstream of the emerging thread. Right Panel: a top view of 
the polymer feedback stress at a distance « H/6 from the bottom wall of the main channel. 


already stressed before, the effect of MV is more pronounced 
with respect to the effect of DV, a conclusion that still holds for 
the Ca and flow parameters used in both the dripping and jet¬ 
ting regime. We therefore choose to report on the effects of vis¬ 
coelasticity in the transition from squeezing to dripping/jetting 
by reporting data only for the case of MV. 

In figure [lO] we report the analysis for the pressure upstream 
of the emerging thread for two different Capillary numbers: 
while for the smaller Capillary number the pressure build-up 
is clearly influenced by viscoelasticity (see also figure [7]), by 
increasing the Capillary number, this effect is less pronounced. 
We remark that the shear viscosity is kept the same in all the 
simulations, so one actually expects to find a less pronounced 
effect of viscoelasticity at larger Capillary numbers, where the 
viscous shear forces start to influence the droplet break-up pro¬ 
cess (see also figure [2]). This is also quantitatively supported 
by the results of Panel (a) of figure ED where we report the 


dimensionless droplet volume V /H 3 as a function of Ca for 
the same values of L 2 considered in the previous figures. The 
flow-rate ratio is kept fixed to Q = 1.0 and Ca is changed in 
the range Ca = 0.001 — 0.03. The Deborah number is ranging 
in the interval De = 2.85 — 5.71 (t p = 2000 — 4000 lbu). At 
increasing the Capillary number, we observe that the tendency 
of viscoelastic stresses to promote a smaller volume soon after 
break-up is somehow less evident. This is also complemented 
by the results in Panel (b) of figure |TT] where we report the 
break-up time T& normalized to the break-up time of the corre¬ 
sponding Newtonian case r£ ewt . Other non trivial effects, how¬ 
ever, are present in the morphology of break-up: while for small 
De we observe that the detachment point shifts downstream of 
the junction (Panel (a) of figure [12]), the increase of the Debo¬ 
rah number favors a stabilization of the break-up point closer to 
the junction (Panel (b) of figure[l2]). Another interesting feature 
found is that viscoelasticity favors the necking process to take 
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Fig. 9: Panels (a)-(b): x and y component of the effective force F eS (see Eq. (V2) ) for a matrix viscoelasticity (MV) case at 
t = ^o + 3.4T shear , De = 3.1, L 2 = 5 x 10 3 , Ca = 0.0026, X = 1.0 and Q = 1.0. We have used the characteristic shear time T shear = H/v c 
as a unit of time, while to is a reference time (the same for all simulations). 


place closer to the channels walls (see Panel (b) in figure [T2]). 
To go deeper into this point, similarly to what we have done in 
figure [9] for the squeezing regime, in figure [13] we analyze the 
effective force in the xy-plane at z = L z /2 for a case with 
L 2 = 5 x 10 3 and Ca = 0.013. Again, an “elastic” region up¬ 
stream of the emerging thread is observed. The straining of the 
fluid upstream of the emerging thread causes a storing of elastic 
energy which is released with an elastic expansion downstream 
of the emerging thread (negative y component of the effective 
force in Panel (b) of figure [13]). This release of elastic energy 
forces the necking process towards the boundary. 

We notice that in the plot of the normalized break-up time 
(Panel (b) of figure ED. a Capillary number of the order of 
Ca = Ca cr « 10 -2 exists, above which the normalized break-up 
time is very close to unity and does not sensibly change with 
De and/or L 2 . We attribute this behaviour to the emergence of 
the jetting regime. The corresponding Newtonian dynamics for 
such Capillary numbers (see figure [2]) indeed reveals that the 
dripping regime is not stable and the droplet detachment point 
gradually moves downstream, until a jet is formed Mm. Sim¬ 
ilarly to figure[l2] density contours of the dispersed phase over¬ 
laid on the polymer feedback stresses at these larger Capillary 
numbers are reported in figure [M] Panel (a) of figure [T4|reports 
the liquid thread just before break-up for a slightly viscoelastic 
case, corresponding to De = 1.42. The break-up point actually 
detaches from the wall as it moves progressively downstream 
(see also figure [2]). Panel (b) of figure [14] reports a case with 
increased Deborah number De = 7.14: we observe that due to 
the presence of the feedback stresses, the break-up point shows 
a slight tendency to move towards the wall, which echoes the 
effects already found in the dripping regime. Notice that due 
to the increase of the Capillary number (Ca = 0.026), the feed¬ 
back stresses are more intense than situations at smaller Ca. 
The small effects on droplet size observed at the higher Capil¬ 
lary numbers in our simulations somehow echo the numerical 


work by Shonibare et al. j55l on T-junctions with viscoelastic 
phases. In particular, Shonibare et al. used 2d numerical simu¬ 
lations, using the Volume of Fluid (VOF) method, to predict the 
size and detachment point of a viscoelastic droplet in a Newto¬ 
nian Matrix. The authors explored both pressure driven flows 
as well as plane Couette flows in the continuous phase: for the 
Newtonian problem they report smaller droplet sizes when the 
cross shear rate is increased, also in agreement with experimen¬ 
tal work m. However, the introduction of viscoelasticity was 
found to have minimal effects on the droplet size. In compar¬ 
ison to our work, some issues are worth being mentioned and 
discussed. Our numerical simulations on droplet viscoelasticity 
(data only partially shown, see also section HD acknowledge 
a small effect on the droplet size as well. As already stressed 
earlier, we find that the elastic effects are more pronounced 
with matrix viscoelasticity and sensibly perturb the droplet size 
when the dispersed phase obstructs the main channel. Even if 
we were dealing with matrix viscoelasticity in the geometry 
of Shonibare et al. isa, we believe that the elastic effects that 
we discussed in section 3.1 would be sensibly reduced as well, 
since the geometry of Shonibare et al. ll55l does not allow a 
considerable obstruction of the main channel, but rather trigger 
droplet detachment and pinch-off based on forces generated by 
the cross-shear rate. 

Another point to be discussed is the importance of the dimen¬ 
sionality in our numerical simulations. In need of an extensive 
study to quantify the importance of the various model parame¬ 
ters, we preliminarly explored the possibility to use two dimen¬ 
sional numerical simulations. In some situations, when the vis¬ 
cosity ratio X is smaller than one and moderate flow rate ratios 
are considered, 2d break-up is actually found to quantitatively 
well compare with 3d break-up (see also figure [3]). However, 
other 2d numerical simulations with viscosity ratio of order 1 
did not capture the underlying physics quantitatively: when the 
squeezing process is about to conclude, long filaments may be 
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stabilized at the detachment point where the side channel meets 
the main channel, thus producing a break-up dynamics which 
is quantitatively different in 2d and 3d. This is a pathology of 
the 2d model, possibly related to the stability of filaments in 2d 
which would be absent in a 3d simulation. Other studies ini 
based on LBM in 2d do not report on droplet break-up with 
those parameters where we observe such pathology. These facts 
said, and not to spoil the correctness of the 3d case, we decided 
to carry out simulations in 3d, while leaving a detailed com¬ 
parison between the 2d and 3d simulations to a future study, 
possibly to identify the correct range of parameters where both 
can be matched 


4 Conclusions 

Microfluidic technologies offer the possibility to generate 
small fluid volumes of dispersed phases (droplets) in con¬ 
tinuous phases. One of the most common droplet generator 
is represented by T-junction geometries mm, where the 
dispersed phase is injected perpendicularly into the main 
channel and the break-up process is induced by forces created 
by the cross-flowing continuous phase. The confinement that 
naturally accompanies these devices has an impact on droplet 
deformation and break-up, which are significantly different 
from those of unbounded droplets. The situation is further 
complicated by the complex properties of the bulk phases, 
whenever constituents have a viscoelastic - rather than New¬ 
tonian - nature. In this paper we have investigated the effects 
of viscoelasticity on the dynamics and break-up of droplets 
in microfluidic T-junctions using numerical simulations of 
dilute polymer solutions at small Capillary numbers up to 
Ca^ 3 x 10 2 and moderate flow-rate ratios Q ~ 1). Our 

numerical model builds upon our previous studies H421I431I and 
is based on a Navier-Stokes (NS) description of the solvent 
based on the lattice Boltzmann models (LBM) coupled to 
constitutive equations for finite extensible non-linear elastic 
dumbbells with the closure proposed by Peterlin (FENE-P 
model). We have used three-dimensional simulations to 
characterize the various characteristic mechanisms of breakup 
in the confined T-junction. Moreover, the various model 
parameters of the FENE-P constitutive equations, including 
the polymer relaxation time t p and the finite extensibility 
parameter L 2 , have been changed to provide quantitative 
details on how the dynamics and break-up properties are 
affected by viscoelasticity, in cases where the viscoelastic 
properties are confined in the dispersed (d) phase (Droplet 
Viscoelasticity, DV), as well as cases where the viscoelastic 
properties are confined in the continuous (c) phase (Matrix 
Viscoelasticity, MV). At fixed flow conditions (i.e. the same 
Ca and Q ) we find that the effects of viscoelasticity are more 
pronounced in the case with MV, which is quantitatively 
attributed to the fact that the flow driving the break-up process 
upstream of the emerging thread can be sensibly perturbed by 
the polymer feedback stresses. This has been evidenced by the 
analysis of simultaneous view of the droplet shape just before 
break-up and the polymer feedback stresses that develop 
in the non-Newtonian phase. In particular, the numerical 
simulations are crucial to elucidate the relative importance of 
the free parameters in the FENE-P model, and to visualize the 


distribution of the polymer feedback stresses. Thanks to these 
insights, it was possible to correlate the distribution of the 
stresses to the corresponding break-up morphology. We also 
tried some preliminary numerical simulations with viscoelastic 
behaviour simultaneously present in both the continuous and 
dispersed phase, and these seem to produce similar effects as 
the matrix viscoelasticity case, at least for the geometry and 
parameters that we explored in our simulations. 

For future investigations, it is surely warranted a comple¬ 
mentary study to highlight the role of viscoelasticity on the 
break-up properties of confined droplets in flow-focusing 
geometries mm. Complementing the available ex¬ 
perimental results with the help of numerical simulations 
would be of extreme interest. Simulations can indeed be used 
to perform in-silico comparative studies, at changing the 
model parameters, to shed lights on the complex properties of 
viscoelastic flows in confined geometries. 
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A Hybrid Lattice Boltzmann Models (LBM) - 
Finite Difference Scheme for dilute Polymer 
solutions 

In this appendix we report the essential details of the numerical 
scheme used. We refer the interested reader to a dedicated pa¬ 
per ED where more extensive technical details are reported, 
together with benchmarks on the rheology of dilute homo¬ 
geneous solutions (including steady shear flow, elongational 
flows, transient shear and oscillatory flows) and viscoelastic 
droplet deformation in confined geometries. We use a hybrid 
algorithm combining a multicomponent Lattice-Boltzmann 
model (LBM) with Finite Differences (FD) schemes, the for¬ 
mer used to model the macroscopic hydrodynamic equations, 
and the latter used to model the polymer dynamics. The LBM 
equations evolve in time the discretized probability density 
function f a i{r,t) to find at position r and time t a fluid par¬ 
ticle of component a = A, B with velocity c;. The dispersed (d) 
and the continuous (c) Newtonian phases in Eqs. Q and ([3]) 
are characterized by a majority of one of the two components, 
i.e. majority of A ( B ) in the dispersed (continuous) phase. The 
LBM evolution scheme with a unitary time-step reads as fol¬ 
lows (56): 

fai{r + C U t + \)-f ai {r,t) = ^ij{faj-faf)+Ki- ( 13 ) 

j 

The collisional operator in the rhs of Eq. (]~3] ) is linear and ex¬ 
presses the relaxation of f a i towards the local equilibrium . 
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(a) Ca = 0.0026 


(b) Ca = 0.0052 


Fig. 10: Analysis of the pressure (P) versus time for different Ca. The pressure is computed upstream of the T-junction and 
Po is a constant reference pressure computed in the static case. Panel (a): we report the normalized pressure versus time with 
fixed Ca = 0.0026 for the Newtonian case (black squares) and two cases with matrix viscoelasticity (MV) at fixed De = 5.7 and 
different L 2 : L 2 = 100 (red circles) and L 2 = 1000 (blue triangles). Panel (b): same as panel (a) for Ca = 0.0052. In all cases we 
have used the characteristic shear time T shear = H/v c as a unit of time, while t ohs is the time when the thread starts to obstruct the 
channel. 
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(a) Matrix Viscoelasticity (MV), Droplet Size 



Fig. 11: Quantitative analysis of the break-up process for different Ca. Panel (a): we report the dimensionless droplet volume 
V/H 3 soon after break-up for a case with matrix viscoelasticity (MV). We choose the polymer relaxation time and the finite 
extensibility parameter L 2 ranging in the interval T p = 2000 — 4000 lbu (De = 2.85 — 5.71 based on (U) and L 2 = 10 2 - 10 3 , 
respectively. The flow rate ratio is kept fixed to Q = 1.0 and the Capillary number is changed in the range Ca = 0.001 — 0.03. 
Panel (b): we report the break-up time T 1> normalized to the break-up time of the corresponding Newtonian case r/ ewt . 


We use the D3Q19 model with 19 velocities 


C/ — < 


(0,0,0) i = 0 

(±1,0,0),(0,±1,0),(0,0,±1) i= 1...6 

(±1, ±1,0), (±1,0,±1), (0, ±1, ±1) f = 7... 18 


(14) 

The expression for the equilibrium distribution is a result of the 
projection onto orthogonal polynomials 15711581 


Jai ~ w iPa 


1 


UCi 


uu : (ciCi-t) 


2 cf 


(15) 


and the weights w; are 



1/3 

1/18 

1/36 


i = 0 

i= 1... 6 , 
z = 7... 18 


(16) 


where c s is the isothermal speed of sound (a constant in the 
model) and u is the fluid velocity. The operator 2^ in Eq. f]~3] ) 
is the same for both components and is characterized by a di¬ 
agonal representation in the mode space : the basis vectors H^ 
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(a) t = to + 11.3T S hear, De = 1.42, L 2 = 10 3 (b) t = to + 10.9T S hear, De = 7.14, L 2 = 10 3 


Fig. 12: Panels (a)-(b): density contours of dispersed phase overlaid on the polymer feedback stress in the stream-flow direction 
(see Eqs. ([!]) and (TT\ ) for a case with matrix viscoelasticity (MV) with L 2 = 10 3 and two different values of De: De = 1.42 
(Panel (a)) and De = 7.14 (Panel (b)). The other parameters are kept fixed to Ca = 0.013, A = 1.0, Q = 1.0. In all cases we have 
used the characteristic shear time T shear = H/v c as a unit of time, while to is a reference time (the same for all simulations). Notice 
that the colorbar of the feedback stress CD is the same. 
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(a) t = to + 8.8 T sh ear; De - 7.14, L 2 = 10 3 (b) t=to+ 8.8 Tshear? De = 7.14, L 2 = 10 3 


Fig. 13: Panels (a)-(b): x and y component of effective force F eff (see Eq. (Y2\ ) for a matrix viscoelasticity (MV) case at t = 
to + 8.8T shear , De = 7.14, L 2 = 10 3 , Ca = 0.013, A = 1.0 and Q = 1.0. At the break-up point we see that there is a net effective 
force which forces the necking process towards the boundary. We have used the characteristic shear time T shear = H/v c as a unit 
of time, while to is a reference time (the same for all simulations). 



(a) t = to + 6.8T s hear; De = 1.42, L 2 = 10 3 (b) t = to + 6.8T s hear; De = 7.14, L 2 = 10 3 


Fig. 14: Panels (a)-(b): density contours of dispersed phase overlaid on the polymer feedback stress in the stream-flow direction 
(see Eqs. 0 and for a case with matrix viscoelasticity (MV) with L 2 = 10 3 and two different values of De: De = 1.42 
(Panel (a)) and De = 7.14 (Panel (b)). The other parameters are kept fixed to Ca = 0.026, A = 1.0, Q = 1.0. In all cases we have 
used the characteristic shear time T shear = H/v c as a unit of time, while to is a reference time (the same for all simulations). Notice 
that the colorbar of the feedback stress GD is the same. 


(k = 0,..., 18) of such mode space are constructed by orthogo- 
nalizing polynomials of the dimensionless velocity vectors [57. 

The basis vectors are used to calculate a complete 
set of moments, the so-called “modes” m a k = Y,i H^ifai (k = 
0,..., 18). The lowest order modes are related to the hydrody¬ 
namic variables, in particular the density (of both components 
and the total one), p a = m a o = Lifai, P = La "'«o = LaPa , 
while the next three moments rh a = {m a \,m a 2 ,m a 0 , are re¬ 


lated to the velocity of the mixture 



m 0 


2 P 


r a i 


2p 


(17) 


The higher order modes refer to the shear and bulk modes in 
the viscous stress tensor, and also other modes (the so called 
“ghost modes”) which do not appear at the level of hydrody¬ 
namic equations. The operator possesses a diagonal repre¬ 
sentation in mode space, hence the collisional term describes 
a linear relaxation of the modes, m p °^ — (1 + A k) m ak + m8 ak > 
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where the “post” indicates the post-collisional mode and where 
the relaxation frequencies —A^ are related to the transport co¬ 


efficients of the modes. The term 


l ak 


is the k-th moment of 


the forcing source A* ■ which embeds the effects of a forcing 
term with density g a H57II59I . The term g = g a in Eq. (T7] ) 
is the total (internal+extemal) force. Forces transfer an amount 
g a of total momentum to the fluid in one time step. The forc¬ 
ing term is determined in such a way that the hydrodynamic 
Eqs. ([20])-(|2T]) are obtained lf6Tl 


Ag : [ciCi-c]!) , 

( 18 ) 

G= 2 + Aj ^ ug + ( ug y + '±+h-t{ u .g) 

(19) 

where the relaxation frequencies of the momentum (—Am), 
bulk (—A^) and shear (—A?) modes appear. LBM is able to re¬ 
produce the continuity equations and the NS equations for the 
total momentum I57ll59ll60l 


as -W( 2 + Am\ 
c]{ 2 


9a 


Wi 

’ c i + -n 

Cf 


dtPa + V-(p a u)=V-D a , (20) 


p [d t u + (u • V)m] = — Vp 


r/ s ( Vu + (Vu) T - -l(V-u) J +rjfol(V-u) 


+ 9 

( 21 ) 


where we have indicated with r/ s , V)^ the shear and bulk viscosi¬ 
ties, respectively. In Eq. ([21]), p = p a = c 2 s p a represents 
the internal (ideal) pressure of the mixture. The quantity D a 
represents the inter-diffusion flux 


Da = 9 




( 22 ) 


with /I a mobility parameter. As for the internal forces, we 
will use the “Shan-Chen” model i44ll62ll63l for multicompo¬ 
nent fluids 


9a(r) =-gABpa(r)Y, E w ‘Pa'( r + c i) c i a,a'=A,B 

i a'^a 

(23) 

where gAB is a parameter that regulates the interactions be¬ 
tween the two components. When gAB is sufficiently large, the 
model can describe stable interfaces with a positive surface ten¬ 
sion. The effect of interaction forces is to introduce an “inter¬ 
action” pressure tensor pM EH, which modifies the internal 
pressure, i.e. P = pt + 

p( M \r) = ^gABpA(r)Y,WiPB(r + Ci)cjCj 

, ' (24) 

+ ~gABPB (f ) E WjpA (r + Ci)c,-Cj. 

Z i 

A tuning of the density in contact with the wall allows for the 
modelling of the wetting properties H46U47II . 

With regard to the transport coefficients of hydrodynamics, the 
relaxation frequencies of the momentum is related to the mo¬ 
bility coefficient 



while the relaxation frequencies of the bulk and shear modes in 
(13\ are related to the viscosity coefficients as 


Vs = ~pc1 


1 1 

T + 2 


j?fc ““3 pc ui + ^ 


(26) 


The numerical simulations presented feature gAB = 1.5 lbu in 
( [23] ), corresponding to a surface tension <7 = 0.1 lbu and associ¬ 
ated bulk densities pA = 2.0 lbu and Pb = 0.1 lbu in the A-rich 
phase. The relaxation frequencies in (26) are set to Am = —1.0 
lbu and A* = A&, thus reproducing the viscous stress tensor 
given in Eqs. 0 and ([3])). The viscosity ratio of the LBM fluid 
is changed by allowing A* to depend on space 


-pc s 


Y + \) =r ls= f?rf(/+(0)) + J?c(/-(0)) 


(27) 


where (/) = (j>(r) 
sen as 


= (p^(r)+pB(r)) - The functions /±(0) are cho- 



l±tanh (0/g) ^ 


(28) 


which allows to recover the Newtonian part of the NS equa¬ 
tions reported in Eqs. 0 and ([3]) with shear viscosities r y and 
Vj c . The smoothing parameter q is chosen sufficiently small so 
as to match analytical predictions on droplet deformation in 
presence of viscoelastic stresses (see |[43l for all details). 

As for the polymer evolution given in Eq. 0, we follow the 
two References 16511661 to solve the FENE-P equation. The 
polymer stress f{r P yg is computed from the FENE-P evo¬ 
lution equation and used to change the shear modes of the 
LBM 1431l57ll58l l. The feedback of the polymers is modu¬ 
lated |[36l in space with the functions f±(<j>) 


p [d t u+ {u- V)tx] = —VP 

+ V [(VdM<l>) + VcfAmVu + (Vrr) r )] (29) 

+ ^V[/M^/ ± (0)]. 

T p 

By using /_ (0), we recover a case where the viscoelastic prop¬ 
erties are confined in the continuous (c) phase, while the use 
of the function /+(0) allows to recover a case where the vis¬ 
coelastic properties are confined in the dispersed (d) phase. 
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